copilot : NeurIPS 2023 workshop :¶

In [47]:
#bbox_of_interest = [-59.958, -9.818, -58.633 , -8.946 ]
#search = catalog.search(collections=["io-biodiversity"], bbox=bbox_of_interest)
#item = next(search.get_items())
#item

everest = [86.9250, 27.9881]
seattle = [-122.332, 47.6062]
grand_canyon = [-112.107676, 36.101690]
mount_fuji = [138.7274, 35.3606]
mont_blanc = [6.865000, 45.832778]
imjalake = [ 86.92586485793316, 27.899354509705653]








areas_of_interestAlaska = {
    "type": "Polygon",
    "coordinates": [
        [
            [-148.56536865234375, 60.80072385643073],
            [-147.44338989257812, 60.80072385643073],
            [-147.44338989257812, 61.18363894915102],
            [-148.56536865234375, 61.18363894915102],
            [-148.56536865234375, 60.80072385643073],
        ]
    ],
}

areas_of_interestImjalake = {
    "type": "Polygon",
    "coordinates": [
        [
        [
              86.92228070099549,
              27.9096329170643
            ],
            [
              86.92228070099549,
              27.882026779073712
            ],
            [
              86.95880774284717,
              27.882026779073712
            ],
            [
              86.95880774284717,
              27.9096329170643
            ],
            [
              86.92228070099549,
              27.9096329170643
            ]
        ]
    ],
}


areas_of_interestLhonaklake = {
    "type": "Polygon",
    "coordinates": [
        [
             [
              88.31558176501454,
              27.97312474305086
            ],
            [
              88.31558176501454,
              27.937831970683348
            ],
            [
              88.34995817854201,
              27.937831970683348
            ],
            [
              88.34995817854201,
              27.97312474305086
            ],
            [
              88.31558176501454,
              27.97312474305086
            ]
        ]
    ],
}


areas_of_interestLhonaklaketoDam = {
    "type": "Polygon",
    "coordinates": [
        [
            [
              88.31217526255193,
              27.960400536405857
            ],
            [
              88.31217526255193,
              27.744973669773245
            ],
            [
              88.57021256641684,
              27.744973669773245
            ],
            [
              88.57021256641684,
              27.960400536405857
            ],
            [
              88.31217526255193,
              27.960400536405857
            ]
        ]
    ],
}





areas_of_interestImjalakeold = {
    "type": "Polygon",
    "coordinates": [
         [

           [
              86.90840990978438,
              27.914442874882738
            ],
            [
              86.90840990978438,
              27.87583921087139
            ],
            [
              86.96126773837943,
              27.87583921087139
            ],
            [
              86.96126773837943,
              27.914442874882738
            ],
            [
              86.90840990978438,
              27.914442874882738
            ]
         ]
    ],
}


areas_of_interestImjalakeGL = {
    "type": "Polygon",
    "coordinates": [
         [
         [
              88.14439581658235,
              27.929979855917935
            ],
            [
              88.14439581658235,
              27.90157666024274
            ],
            [
              88.21527228069135,
              27.90157666024274
            ],
            [
              88.21527228069135,
              27.929979855917935
            ],
            [
              88.14439581658235,
              27.929979855917935
            ]

         ]
    ],
}


   



#areas_of_interest = {"type": "Point", "coordinates": imjalake}

#areas_of_interest = areas_of_interestImjalake

#areas_of_interest = areas_of_interestLhonaklake

areas_of_interest = areas_of_interestImjalake


time_of_interest = "2019-06-01/2023-10-01"
time_of_interest2 = "2021-08-01/2022-12-01"

https://geojson.io/#map=11.12/27.8872/86.9057

In [48]:
yearofinterest = str(2021)
time_of_interest = yearofinterest + "-07-01/" + yearofinterest + "-08-01"
time_of_interest
Out[48]:
'2021-07-01/2021-08-01'
In [49]:
import rasterio
from rasterio import windows
from rasterio import features
from rasterio import warp

import numpy as np
from PIL import Image

def displayRawImage(least_cloudy_item):
    asset_href = least_cloudy_item.assets["visual"].href

    with rasterio.open(asset_href) as ds:
        aoi_bounds = features.bounds(areas_of_interest)
        print("area of interest")
        print(aoi_bounds)
        warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
        print("wraped aoi bounds")
       # print(warped_aoi_bounds)    
        aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
        print(aoi_window)    
        band_data = ds.read(window=aoi_window)
        #print(band_data)

    img = Image.fromarray(np.transpose(band_data, axes=[1, 2, 0]))
    w = img.size[0]
    h = img.size[1]
    aspect = w / h
    target_w = 800
    target_h = (int)(target_w / aspect)
    return img.resize((target_w, target_h), Image.Resampling.BILINEAR)
In [50]:
print(time_of_interest2)

2021-08-01/2022-12-01
In [51]:
def displayRawImagesentinel1(least_cloudy_item):

    img_arr = rioxarray.open_rasterio(least_cloudy_item.assets["rendered_preview"].href).squeeze()
    #print(item.datetime.strftime("%x"))
    img_arr.plot.imshow(figsize=(25, 20));
    return
In [52]:
from IPython.display import Image
In [54]:
from pystac.extensions.eo import EOExtension as eo
import pystac_client
import planetary_computer

import pystac_client
import planetary_computer

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1/",
    modifier=planetary_computer.sign_inplace,
)
In [55]:
import numpy as np
from PIL import Image

import rioxarray
import matplotlib.pyplot as plt
In [56]:
def displayRawImage222222(least_cloudy_item):
    asset_href = least_cloudy_item.assets["visual"].href

    with rasterio.open(asset_href) as ds:
        aoi_bounds = features.bounds(areas_of_interest)
        print("area of interest")
        print(aoi_bounds)
        warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
        print("wraped aoi bounds")
       # print(warped_aoi_bounds)    
        aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
        print(aoi_window)    
        band_data = ds.read(window=aoi_window)
        #print(band_data)
    return band_data[2]
In [ ]:
In [ ]:
In [57]:
def displayRawImageSAC(least_cloudy_item):
    asset_href = least_cloudy_item.assets["SCL"].href

    with rasterio.open(asset_href) as ds:
        aoi_bounds = features.bounds(areas_of_interest)
        print("area of interest")
        print(aoi_bounds)
        warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
        print("wraped aoi bounds")
       # print(warped_aoi_bounds)    
        aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
        print(aoi_window)    
        band_data = ds.read(window=aoi_window)
        #print(band_data)

        band_datashow = band_data == 6
       # plt.imshow(band_datashow[0])
    return band_data[0]
In [ ]:
In [59]:
for i in range(2012, 2024, 1):
    print(i)
    yearofinterest = str(i)
    time_of_interest = yearofinterest + "-08-01/" + yearofinterest + "-11-01"    
    print(time_of_interest)
    search = catalog.search(collections=["sentinel-2-l2a"], intersects=areas_of_interest, datetime=time_of_interest,   query={"eo:cloud_cover": {"lt": 50, "gt":2 }})
   # search = catalog.search(collections=["sentinel-2-l2a"], intersects=areas_of_interest, datetime=time_of_interest )
    items = search.item_collection()
    print(f"Returned {len(items)} Items")
    
    if (len(items) >0):
        least_cloudy_item = min(items, key=lambda item: eo.ext(item).cloud_cover)
        item = least_cloudy_item
        answer = displayRawImage(least_cloudy_item)
        plt.imshow(answer)
        plt.title(item.datetime.strftime("%B %d, %Y"))
        plt.show()
        answer = displayRawImageSAC(least_cloudy_item)
        print(item.datetime.strftime("%x"))    
        plt.imshow(answer)
        plt.title(item.datetime.strftime("%B %d, %Y"))
      
        plt.show()
        

2012
2012-08-01/2012-11-01
Returned 0 Items
2013
2013-08-01/2013-11-01
Returned 0 Items
2014
2014-08-01/2014-11-01
Returned 0 Items
2015
2015-08-01/2015-11-01
Returned 0 Items
2016
2016-08-01/2016-11-01
Returned 2 Items
area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)

area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496)
10/30/16

2017
2017-08-01/2017-11-01
Returned 3 Items
area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)

area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496)
10/30/17

2018
2018-08-01/2018-11-01
Returned 8 Items
area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)

area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496)
10/30/18

2019
2019-08-01/2019-11-01
Returned 6 Items
area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)

area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496)
10/15/19

2020
2020-08-01/2020-11-01
Returned 7 Items
area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)

area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496)
10/29/20

2021
2021-08-01/2021-11-01
Returned 2 Items
area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)

area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496)
10/14/21

2022
2022-08-01/2022-11-01
Returned 5 Items
area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)

area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496)
10/24/22

2023
2023-08-01/2023-11-01
Returned 2 Items
area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)

area of interest
(86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643)
wraped aoi bounds
Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496)
09/29/23
In [60]:
def displayRawImageSNOW(least_cloudy_item):

    img_arr = rioxarray.open_rasterio(least_cloudy_item.assets["visual"].href).squeeze()
    print(item.datetime.strftime("%x"))
    img_arr.plot.imshow(figsize=(25, 20));
    return
In [61]:
def displayRawImageSNOW2(least_cloudy_item):

    img_arr = rioxarray.open_rasterio(least_cloudy_item.assets["visual"].href).squeeze()
    print(item.datetime.strftime("%x"))
    
    return img_arr
In [63]:
def displayRawImageSNOWCLASSIFY(least_cloudy_item):

    img_arr = rioxarray.open_rasterio(least_cloudy_item.assets["SCL"].href).squeeze()
  #  print(least_cloudy_item.datetime.strftime("%x"))
    #https://custom-scripts.sentinel-hub.com/custom-scripts/sentinel-2/scene-classification/
    img_arrnew = img_arr == 11    #11-snow, 6-water 
    return img_arrnew
In [45]:
def displayRawImageSNOWCLASSIFY2(least_cloudy_item):

    img_arr = rioxarray.open_rasterio(least_cloudy_item.assets["SCL"].href).squeeze()
  #  print(least_cloudy_item.datetime.strftime("%x"))
    #https://custom-scripts.sentinel-hub.com/custom-scripts/sentinel-2/scene-classification/
    img_arrnew = img_arr == 11    #11-snow, 6-water 
    return img_arr
In [23]:
from datashader.transfer_functions import shade, stack
from datashader.colors import Elevation
from xrspatial import hillshade

import pystac_client
import planetary_computer
import rioxarray


search = catalog.search(collections=["nasadem"], intersects=areas_of_interest)

items = search.item_collection()

print(f"Returned {len(items)} Items")

item = items[0]

da = (
    rioxarray.open_rasterio(item.assets["elevation"].href)
    .squeeze()
    .drop("band")[:-1, :-1]
    .coarsen({"y": 5, "x": 5})
    .mean()
)

Returned 1 Items
In [24]:
from datashader.transfer_functions import shade
import matplotlib.pyplot as plt
from xrspatial.classify import natural_breaks


for i in range(1):
    print(i)
    yearofinterest = str(i)
    time_of_interest = yearofinterest + "-06-01/" + yearofinterest + "-08-01"    
    print(time_of_interest)
    #search = catalog.search(collections=["nasadem"], intersects=areas_of_interest,    query={"eo:cloud_cover": {"lt": 450, "gt":2 }})
    search = catalog.search(collections=["nasadem"], intersects=areas_of_interest  )
    items = search.item_collection()
    print(f"Returned {len(items)} Items")
    
    if (len(items) >0):
        item = items[0]
        answer = rioxarray.open_rasterio(item.assets["elevation"].href).squeeze()
        
        natural_breaks_result = natural_breaks(answer, num_sample=20000, k=15)
        shade(natural_breaks_result, cmap=plt.get_cmap("terrain"), how="linear")
    
 
        
# classify the image using natural breaks

0
0-06-01/0-08-01
Returned 1 Items
In [ ]:
least_cloudy_item = min(items, key=lambda item: eo.ext(item).cloud_cover)

print(
    f"Choosing {least_cloudy_item.id} from {least_cloudy_item.datetime.date()}"
    f" with {eo.ext(least_cloudy_item).cloud_cover}% cloud cover"
)
In [ ]:
item = least_cloudy_item
In [ ]:
In [ ]:
asset_href = least_cloudy_item.assets["visual"].href
In [ ]:
items = list(search.items())
for item in items:
    print(item)
In [ ]:
import rich.table

asset_table = rich.table.Table("Asset Key", "Asset Title")
for key, value in items[-1].assets.items():
    asset_table.add_row(key, value.title)
asset_table
In [35]:
import rasterio
from rasterio import windows
from rasterio import features
from rasterio import warp

import numpy as np
from PIL import Image

asset_href = least_cloudy_item.assets["visual"].href

with rasterio.open(asset_href) as ds:
    aoi_bounds = features.bounds(areas_of_interest)
    print("area of interest")
    print(aoi_bounds)
    warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
    print("wraped aoi bounds")
    print(warped_aoi_bounds)    
    aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
    print(aoi_window)    
    band_data = ds.read(window=aoi_window)
    print(band_data)
   

area of interest
(88.14439581658235, 27.90157666024274, 88.21527228069135, 27.929979855917935)
wraped aoi bounds
(612596.7648950818, 3086825.9764255327, 619602.1054647851, 3090039.936055328)
Window(col_off=1259.6764895081797, row_off=998.006394467142, width=700.534056970333, height=321.39596297958633)
[[[198 255 255 ... 182 193 200]
  [212 255 255 ... 189 202 202]
  [245 255 255 ... 197 208 206]
  ...
  [255 255 255 ...  64  69  59]
  [255 255 255 ...  75  72  60]
  [255 255 255 ...  80  75  69]]

 [[204 255 255 ... 161 168 171]
  [182 255 255 ... 166 177 176]
  [171 255 255 ... 174 184 181]
  ...
  [255 255 255 ...  69  67  58]
  [255 255 255 ...  83  69  64]
  [255 255 255 ...  82  80  94]]

 [[166 255 255 ... 117 124 125]
  [141 255 255 ... 120 129 127]
  [113 255 255 ... 126 135 132]
  ...
  [255 255 255 ...  73  65  58]
  [255 255 255 ...  77  78  87]
  [255 255 255 ...  74  86  94]]]
In [36]:
areas_of_interest
Out[36]:
{'type': 'Polygon',
 'coordinates': [[[88.14439581658235, 27.929979855917935],
   [88.14439581658235, 27.90157666024274],
   [88.21527228069135, 27.90157666024274],
   [88.21527228069135, 27.929979855917935],
   [88.14439581658235, 27.929979855917935]]]}
In [ ]:
In [38]:
import rasterio
from rasterio import windows
from rasterio import features
from rasterio import warp

import numpy as np
from PIL import Image

asset_href = least_cloudy_item.assets["SCL"].href

with rasterio.open(asset_href) as ds:
    aoi_bounds = features.bounds(areas_of_interest)
    print(aoi_bounds)
    warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
    print(warped_aoi_bounds)    
    aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
    print(aoi_window)    
    band_data = ds.read(window=aoi_window)
    print(band_data)
   

(88.14439581658235, 27.90157666024274, 88.21527228069135, 27.929979855917935)
(612596.7648950818, 3086825.9764255327, 619602.1054647851, 3090039.936055328)
Window(col_off=629.8382447540898, row_off=499.003197233571, width=350.2670284851665, height=160.69798148979316)
[[[11 11 11 ...  5  5  5]
  [11 11 11 ...  5  5  5]
  [11 11 11 ...  5  5  5]
  ...
  [11 11 11 ...  7  2  2]
  [11 11 11 ...  7  7  2]
  [11 11 11 ...  7  7  7]]]
In [39]:
 band_data.shape
Out[39]:
(1, 161, 350)
In [40]:
type(band_data)
Out[40]:
numpy.ndarray

This item's data can be read in using xarray's open_rasterio function to load the data into an array:

In [43]:
import rioxarray

item = least_cloudy_item

img_arr = rioxarray.open_rasterio(item.assets["visual"].href).squeeze()
print(item.datetime.strftime("%x"))
img_arr.plot.imshow(figsize=(25, 20));

08/12/22
In [46]:
img_arr
Out[46]:
<xarray.DataArray (y: 5490, x: 5490)>
array([[11, 11, 11, ...,  5,  5,  5],
       [11, 11, 11, ...,  5,  5,  5],
       [11, 11, 11, ...,  5,  5,  5],
       ...,
       [ 9,  9,  9, ...,  4,  4,  4],
       [ 9,  9,  9, ...,  4,  4,  4],
       [ 9,  9,  9, ...,  4,  4,  4]], dtype=uint8)
Coordinates:
    band         int64 1
  * x            (x) float64 6e+05 6e+05 6e+05 ... 7.098e+05 7.098e+05 7.098e+05
  * y            (y) float64 3.1e+06 3.1e+06 3.1e+06 ... 2.99e+06 2.99e+06
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    _FillValue:     0
    scale_factor:   1.0
    add_offset:     0.0
xarray.DataArray
  • y: 5490
  • x: 5490
  • 11 11 11 11 11 11 5 5 5 5 5 11 11 11 ... 4 4 4 4 4 4 4 4 4 4 4 4 4 4
    array([[11, 11, 11, ...,  5,  5,  5],
       [11, 11, 11, ...,  5,  5,  5],
       [11, 11, 11, ...,  5,  5,  5],
       ...,
       [ 9,  9,  9, ...,  4,  4,  4],
       [ 9,  9,  9, ...,  4,  4,  4],
       [ 9,  9,  9, ...,  4,  4,  4]], dtype=uint8)
    • band
      ()
      int64
      1
      array(1)
    • x
      (x)
      float64
      6e+05 6e+05 ... 7.098e+05 7.098e+05
      array([600010., 600030., 600050., ..., 709750., 709770., 709790.])
    • y
      (y)
      float64
      3.1e+06 3.1e+06 ... 2.99e+06
      array([3100010., 3099990., 3099970., ..., 2990270., 2990250., 2990230.])
    • spatial_ref
      ()
      int64
      0
      crs_wkt :
      PROJCS["WGS 84 / UTM zone 45N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32645"]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      horizontal_datum_name :
      World Geodetic System 1984
      projected_crs_name :
      WGS 84 / UTM zone 45N
      grid_mapping_name :
      transverse_mercator
      latitude_of_projection_origin :
      0.0
      longitude_of_central_meridian :
      87.0
      false_easting :
      500000.0
      false_northing :
      0.0
      scale_factor_at_central_meridian :
      0.9996
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 45N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32645"]]
      GeoTransform :
      600000.0 20.0 0.0 3100020.0 0.0 -20.0
      array(0)
    • x
      PandasIndex
      PandasIndex(Index([600010.0, 600030.0, 600050.0, 600070.0, 600090.0, 600110.0, 600130.0,
       600150.0, 600170.0, 600190.0,
       ...
       709610.0, 709630.0, 709650.0, 709670.0, 709690.0, 709710.0, 709730.0,
       709750.0, 709770.0, 709790.0],
      dtype='float64', name='x', length=5490))
    • y
      PandasIndex
      PandasIndex(Index([3100010.0, 3099990.0, 3099970.0, 3099950.0, 3099930.0, 3099910.0,
       3099890.0, 3099870.0, 3099850.0, 3099830.0,
       ...
       2990410.0, 2990390.0, 2990370.0, 2990350.0, 2990330.0, 2990310.0,
       2990290.0, 2990270.0, 2990250.0, 2990230.0],
      dtype='float64', name='y', length=5490))
  • AREA_OR_POINT :
    Area
    _FillValue :
    0
    scale_factor :
    1.0
    add_offset :
    0.0
In [ ]:
In [47]:
img_arrnew = img_arr == 11    #11-snow, 6-water
In [ ]:
In [48]:
img_arrnew
Out[48]:
<xarray.DataArray (y: 5490, x: 5490)>
array([[ True,  True,  True, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]])
Coordinates:
    band         int64 1
  * x            (x) float64 6e+05 6e+05 6e+05 ... 7.098e+05 7.098e+05 7.098e+05
  * y            (y) float64 3.1e+06 3.1e+06 3.1e+06 ... 2.99e+06 2.99e+06
    spatial_ref  int64 0
xarray.DataArray
  • y: 5490
  • x: 5490
  • True True True True True True ... False False False False False False
    array([[ True,  True,  True, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]])
    • band
      ()
      int64
      1
      array(1)
    • x
      (x)
      float64
      6e+05 6e+05 ... 7.098e+05 7.098e+05
      array([600010., 600030., 600050., ..., 709750., 709770., 709790.])
    • y
      (y)
      float64
      3.1e+06 3.1e+06 ... 2.99e+06
      array([3100010., 3099990., 3099970., ..., 2990270., 2990250., 2990230.])
    • spatial_ref
      ()
      int64
      0
      crs_wkt :
      PROJCS["WGS 84 / UTM zone 45N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32645"]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      horizontal_datum_name :
      World Geodetic System 1984
      projected_crs_name :
      WGS 84 / UTM zone 45N
      grid_mapping_name :
      transverse_mercator
      latitude_of_projection_origin :
      0.0
      longitude_of_central_meridian :
      87.0
      false_easting :
      500000.0
      false_northing :
      0.0
      scale_factor_at_central_meridian :
      0.9996
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 45N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32645"]]
      GeoTransform :
      600000.0 20.0 0.0 3100020.0 0.0 -20.0
      array(0)
    • x
      PandasIndex
      PandasIndex(Index([600010.0, 600030.0, 600050.0, 600070.0, 600090.0, 600110.0, 600130.0,
       600150.0, 600170.0, 600190.0,
       ...
       709610.0, 709630.0, 709650.0, 709670.0, 709690.0, 709710.0, 709730.0,
       709750.0, 709770.0, 709790.0],
      dtype='float64', name='x', length=5490))
    • y
      PandasIndex
      PandasIndex(Index([3100010.0, 3099990.0, 3099970.0, 3099950.0, 3099930.0, 3099910.0,
       3099890.0, 3099870.0, 3099850.0, 3099830.0,
       ...
       2990410.0, 2990390.0, 2990370.0, 2990350.0, 2990330.0, 2990310.0,
       2990290.0, 2990270.0, 2990250.0, 2990230.0],
      dtype='float64', name='y', length=5490))
In [49]:
snowcount = sum((sum(img_arrnew)))
print("Snow area = " + str(snowcount))

Snow area = <xarray.DataArray ()>
array(2291118)
Coordinates:
    band         int64 1
    spatial_ref  int64 0
In [50]:
type(snowcount)
Out[50]:
xarray.core.dataarray.DataArray
In [51]:
snowcount.values
Out[51]:
array(2291118)
In [53]:
type(img_arrnew)
Out[53]:
xarray.core.dataarray.DataArray
In [ ]:
In [65]:
import rioxarray

img_arr = rioxarray.open_rasterio(item.assets["elevation"].href).squeeze().drop("band")
img_arr.plot.imshow(figsize=(15, 10));

Y